#include "transient.h"
#define MY_PI 3.14159265359
GHashTable* sources;

void print_func(gpointer key, gpointer g_value, gpointer user_data) {
	printf("(%s)\t", (gchar *)key);
	printf("%s\n", (gchar *)g_value);
}

void TRAN_init()
{
	sources = g_hash_table_new(g_str_hash, g_str_equal);
}

void TRAN_new_source(char *name,char *data)
{
	g_hash_table_insert(sources, g_strdup(name),g_strdup(data));
}

void TRAN_print()
{
	g_hash_table_foreach(sources, print_func, NULL);
}

void TRAN_get_value(char *name, double t, double *ret_val, int k)
{
	const char *dels = " \t()\n=,";
	char *lookup_value = g_hash_table_lookup(sources, g_strdup(name));
	char *tuple = (char *)malloc(strlen(lookup_value));
	strcpy(tuple,lookup_value);
	char *type = (char *)malloc(sizeof(gchar));
	type = strtok(tuple,dels);
	
	char *tok = malloc(sizeof(gchar));

	if(!strcmp(type, "EXP")) {	//(i1 i2 td1 tc1 td2 tc2)
		tok = strtok(NULL,dels);		
		double i1 = strtod(tok,NULL);

		tok = strtok(NULL,dels);
		double i2 = strtod(tok,NULL);
		
		tok = strtok(NULL,dels);
		double td1 = strtod(tok,NULL);
		
		tok = strtok(NULL,dels);
		double tc1 = strtod(tok,NULL);
		
		tok = strtok(NULL,dels);		
		double td2 = strtod(tok,NULL);
		
		tok = strtok(NULL,dels);
		double tc2 = strtod(tok,NULL);
		
		
		if( t <= td1){			
			*ret_val = i1;
		}
		else if(t <= td2){
			*ret_val = i1 + (i2 - i1)*(1 - exp(-(t-td1)/tc1));
		}
		else{
			*ret_val = i1 + (i2 - i1)*(exp( -(t-td2)/tc2) - exp(-(t-td1)/tc1));
		}
	}
	else if(!strcmp(type, "SIN")) {		//i1 ia fr td df ph)
		
		tok = strtok(NULL,dels);		
		double i1 = strtod(tok,NULL);

		tok = strtok(NULL,dels);
		double ia = strtod(tok,NULL);
		
		tok = strtok(NULL,dels);
		double fr = strtod(tok,NULL);
		
		tok = strtok(NULL,dels);
		double td = strtod(tok,NULL);
		
		tok = strtok(NULL,dels);		
		double df = strtod(tok,NULL);
		
		tok = strtok(NULL,dels);
		double ph = strtod(tok,NULL);
		
		if(t <= td){			
			*ret_val = i1 + ia*sin(2*MY_PI*floor(ph/360));
		}
		else{
			*ret_val = i1 + ia*sin(2*MY_PI*fr*(t-td) + 2*MY_PI*floor(ph/360))*exp(-(t-td)/df);
		}
	}
	
	else if(!strcmp(type, "PULSE")) {		//(i1 i2 td tr tf pw per)
		
		tok = strtok(NULL,dels);		
		double i1 = strtod(tok,NULL);

		tok = strtok(NULL,dels);
		double i2 = strtod(tok,NULL);
		
		tok = strtok(NULL,dels);
		double td = strtod(tok,NULL);
		
		tok = strtok(NULL,dels);
		double tr = strtod(tok,NULL);
		
		tok = strtok(NULL,dels);		
		double tf = strtod(tok,NULL);
		
		tok = strtok(NULL,dels);
		double pw = strtod(tok,NULL);

		tok = strtok(NULL,dels);
		double per = strtod(tok,NULL);
	
		if(t <= td){			
			*ret_val = i1;
		}
		else if( (td + (k*per) <= t) && (t <= td + tr +(k*per)) ){
			*ret_val = (i2-i1);
		}
		else if( (td + tr + (k*per) <= t) && (t <= td + tr + pw +(k*per)) ){
			*ret_val = i2;
		}
		else if( (td + tr + pw + (k*per) <= t) && (t <= td + tr + pw + tf + (k*per)) ){
			*ret_val = i1-i2;
		}
		else if( (td + tr + pw + tf + (k*per) <= t) && (t <= td + per + (k*per)) ){
			*ret_val = i1;
		}
	}
	else if(!strcmp(type, "PWL")) {		//i1 ia fr td df ph)
		
		tok = strtok(NULL,dels);		
		double time = strtod(tok,NULL);
		
		tok = strtok(NULL,dels);
		double val = strtod(tok,NULL);
		
		while(time < t)
		{
			tok = strtok(NULL,dels);		
			time = strtod(tok,NULL);
		
			tok = strtok(NULL,dels);
			val = strtod(tok,NULL);
		}
		*ret_val = val;
	}
	else if(!strcmp(type, "DC")){
		
		tok = strtok(NULL,dels);		
		double dc_val = strtod(tok,NULL);	
			
		*ret_val = dc_val;
	}
	//free(tok);
	//printf("(%s)\t%g---%g\n",type,t,*ret_val);	
	free(tuple);
}

void TRAN_backward_euler(const double *mna_g, const double *mna_c, const double *e, const double *x, 
										int n, double h, double *A, double *RHS)
{
	int i,j;
	
	double I[n*n];
	memset(I, 1, n*n*sizeof(double));
	for(i = 0; i < n; ++i)
		I[i*n+i] = 1;
	
	//aristero melos
	double temp[n*n];	
	for(i = 0; i < n; ++i)
		for(j = 0; j < n; ++j)
			A[i*n+j] = mna_g[i*n+j] + 1/h*mna_c[i*n + j];


	//de3i melos
	double temp_b[n];
	cblas_dcopy (n, e, 1, temp_b, 1);
	cblas_dgemv (CblasRowMajor,CblasNoTrans,n,n,1/h,mna_c,n,x,1,1,temp_b,1);
	cblas_dcopy (n, temp_b, 1, RHS, 1);

}

void TRAN_backward_euler_sparse(const cs *mna_g, const cs *mna_c, const double *e, const double *x, 
										int n, double h, cs **A, double *RHS)
{
	//aristero melos
	//cs *res;
	//res = *A;
	*A = cs_add (mna_g, mna_c, 1, 1/h);

	//de3i melos
	double temp_b[n];
	cs *temp;	
	memset(temp_b, 0, n*sizeof(double));
	
	temp = cs_add (mna_g, mna_c, 0, 1/h);	/* C = alpha*A + beta*B */
	cs_gaxpy (temp, x, temp_b);
	
	int i;
	for(i = 0; i < n; i++)
		RHS[i] = e[i] + temp_b[i];

	//cs_print(res,1);

}



//TODO
void TRAN_trapezoidial(const double *mna_g, const double *mna_c, const double *e, const double *e_prev,
						const double *x_prev, int n, double h, double *A, double *RHS)
{
	int i,j;
	double temp_matrix[n*n];
	double I[n*n];

	memset(I, 1, n*n*sizeof(double));
	for(i = 0; i < n; ++i)
		I[i*n+i] = 1;

	
	for(i = 0; i < n; ++i){
		for(j = 0; j < n; ++j){
			A[i*n+j] = mna_g[i*n+j] + 2/h*mna_c[i*n + j];
			temp_matrix[i*n + j] = mna_g[i*n+j] - 2/h*mna_c[i*n + j];
		}
		RHS[i] = e[i] + e_prev[i];
	}
	cblas_dgemv(CblasRowMajor,CblasNoTrans,n,n,-1,	temp_matrix,n,x_prev,1,1,RHS,1);

}
void TRAN_trapezoidial_sparse(const cs *mna_g, const cs *mna_c, const double *e, const double *e_prev,
						const double *x_prev, int n, double h, cs **A, double *RHS)
{
	//aristero melos
	*A = cs_add (mna_g, mna_c, 1, 2/h);

	//de3i melos
	double temp_b[n];
	cs *temp;	
	memset(temp_b, 0, n*sizeof(double));
	
	temp = cs_add (mna_g, mna_c, 1, -2/h);	/* C = alpha*A + beta*B */
	cs_gaxpy (temp, x_prev, temp_b);
	
	int i;
	for(i = 0; i < n; i++)
		RHS[i] = e[i] + e_prev[i] - temp_b[i];
}


